#include "ParticleSource1D_Load.h"
#include "../Parallel/ParallelMPI1D.h"
#include "../Field/Field1D.h"

#include <cmath>
#include <iomanip>
#include <algorithm>
#include <ostream>
#include <sstream>

using namespace std;


//Constructor
ParticleSource1D_Load::ParticleSource1D_Load(PicParams* params_, ParallelMPI* pmpi_, PhysicalField* fields_, vector<Species*>* vecSpecies_, ParticleSourceParameter ps_param_):
    ParticleSource1D (params_, pmpi_, fields_, vecSpecies_, ps_param_)
{
    dt_ov_dx = params->dt / params->cell_length[0];
    dt = params->dt;

    YZArea = 1.0;
    load_temperature = ps_param_.temperature;
    load_temperature_init = ps_param_.temperature;

    if(ps_param_.load_type == "nq")
    {
        is_begin_change_dn = false;
        load_dq = ps_param_.dq;
        //initial load_dn is 5 times larger, to pass the inital phase quickly
        //load_dn = 5.0 * load_dq / (1.5 * load_temperature_init);
        load_dn = ps_param_.load_dn_scale * load_dq / (1.5 * load_temperature_init);
        load_dn_final = load_dq / (1.5 * load_temperature_init);
    }
    else
    {
        load_dn = ps_param_.dn;
        load_dq = ps_param_.dn * ps_param_.temperature;        
    }
    

    //load particles by number_density and Temperature
    if(ps_param_.load_type == "nT")
    {
        for(int i = 0; i < number_particle_in_each_cell.size(); i++)
        {
            number_particle_in_each_cell[i] = ps_param_.number_density / params->species_param[ps_param_.species_number].weight;
            //cout<<"number_particle_in_each_cell "<<number_particle_in_each_cell[ibin]<<endl;
        }
    }
    else if(ps_param_.load_type == "dn")
    {
        if(loadTimeStepVector.size() != 0)
        {
            load_temperature = load_temperatureVector[0];
            load_dn = load_dnVector[0];
        }

        load_step = 1.0 + ps_param_.load_particle_number * params->species_param[ps_param_.species_number].weight / (load_dn * params->dt);
        load_number_double = load_dn * load_step * params->dt / params->species_param[ps_param_.species_number].weight;
        load_number_rem = 0.0;
        load_number_rem_tot = 0.0;
        log_info<<"load_step = "<<load_step;
    }

    else if(ps_param_.load_type == "nq")
    {
        //load_number_init_double to keep dq = ps_param_.dq
        load_step = ps_param_.load_particle_number * params->species_param[ps_param_.species_number].weight / (load_dn * params->dt);
        load_number_init_double = load_dn * load_step * params->dt / params->species_param[ps_param_.species_number].weight;
        load_number_final_double = load_dn_final * load_step * params->dt / params->species_param[ps_param_.species_number].weight;
        load_number_rem_tot = 0.0;
        log_info<<"load_step = "<<load_step;

        int mpiSize = pmpi->getSize();
        if(mpiSize % 2 == 0)
        {
            mpiRank_source_middle = mpiSize / 2;
            index_source_middle = 2;
        }
        else
        {
            mpiRank_source_middle = mpiSize / 2;
            index_source_middle = params->dims_space[0] / 2;
        }
    }


    //the MPI domain is not in the source region
    if(pmpi->getDomainLocalMax(0) <= ps_param_.load_position_start[0] || pmpi->getDomainLocalMin(0) >= ps_param_.load_position_end[0])
    {
        load_cell_start[0] = 0;
        load_cell_end[0] = 0;
    }
    //the left end of source region is in the MPI domain
    else if(pmpi->getDomainLocalMin(0) <= ps_param_.load_position_start[0] && pmpi->getDomainLocalMax(0) > ps_param_.load_position_start[0]
    && pmpi->getDomainLocalMax(0) <= ps_param_.load_position_end[0])
    {
        double temp = (ps_param_.load_position_start[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        load_cell_start[0] = temp + 0.001;
        if(load_cell_start[0] < 0) { load_cell_start[0] = 0; }
        load_cell_end[0] = params->dims_space[0] - 1;
    }
    //the right end of source region is in the MPI domain
    else if(pmpi->getDomainLocalMin(0) < ps_param_.load_position_end[0] && pmpi->getDomainLocalMax(0) >= ps_param_.load_position_end[0]
    && pmpi->getDomainLocalMin(0) >= ps_param_.load_position_start[0])
    {
        load_cell_start[0] = 0;
        load_cell_end[0] = (ps_param_.load_position_end[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        if(load_cell_end[0] > params->dims_space[0] - 1) { load_cell_end[0] = params->dims_space[0] - 1; }
    }
    //the whole MPI domain is in the source region
    else if(pmpi->getDomainLocalMin(0) >= ps_param_.load_position_start[0] && pmpi->getDomainLocalMax(0) <= ps_param_.load_position_end[0])
    {
        load_cell_start[0] = 0;
        load_cell_end[0] = params->dims_space[0] - 1;
    }
    //the whole MPI domain source region is in the bin
    else if(pmpi->getDomainLocalMin(0) < ps_param_.load_position_start[0] && pmpi->getDomainLocalMax(0) > ps_param_.load_position_end[0])
    {
        double temp = (ps_param_.load_position_start[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        load_cell_start[0] = temp + 0.001;
        if(load_cell_start[0] < 0) { load_cell_start[0] = 0; }
        load_cell_end[0] = (ps_param_.load_position_end[0] - pmpi->getDomainLocalMin(0)) / params->cell_length[0];
        if(load_cell_end[0] > params->dims_space[0] - 1) { load_cell_end[0] = params->dims_space[0] - 1; }
    }

    if(load_cell_end[0] != load_cell_start[0])
    {
        log_debug<<"load start and end: "<<load_cell_start[0]<<" "<<load_cell_end[0];
    }

    //Parameters for "nq"
    load_source_density_pre = 0.0;

    //Parameters for "dn"
    nextTimeStep = 0;
    nextTimeStep_index = 1;


}

ParticleSource1D_Load::~ParticleSource1D_Load()
{

}


void ParticleSource1D_Load::emitLoad(int itime)
{
    Species1D   *species1D;

    //Parameters for "nq"
    double zoom_factor;
    double load_dn_temp;

    species1D = static_cast<Species1D*>((*vecSpecies)[ps_param.species_number]);


    //log_info<<<<"particle source 1d load start";
    if(ps_param.is_every_time == false && itime > 1) { return; }

    //for load_type = nT
    if(itime % ps_param.timestep_update != 0){ return; }

    //for load_type = nq: n(number_density) and q(heat power) in source region is fixed, q = n * temperature from input file
    //dn is dynamicly changed

    if(ps_param.load_type == "nq")
    {
        if(itime%load_step == 0)
        {
            //cout<<"nq begin"<<endl;
            Field1D<double>* rho1D = static_cast<Field1D<double>*>(fields->rho_s[ps_param.species_number_dependent]);
            if(pmpi->getRank() == mpiRank_source_middle)
            {
                load_source_density = (*rho1D)(index_source_middle);
            }
            pmpi->bcast_double(&load_source_density, 1, mpiRank_source_middle);


            if(is_begin_change_dn == false && load_source_density > ps_param.dn_scale_down * ps_param.number_density)
            {
                double load_dn_init = load_dn;
                load_dn =  0.5 * ps_param.number_density * sqrt(ps_param.temperature * params->const_e / 3.34524316e-27) / (ps_param.load_position_end[0] - ps_param.load_position_start[0]);
                load_dn_pid = load_dn;
                if(load_dn > load_dn_init)
                {
                    load_dn = load_dn_init;
                }                
                is_begin_change_dn = true;
                load_number_init_double = load_number_final_double;

                i_value = 0.0;
            }

            if(is_begin_change_dn == true)
            {
                /*
                if(load_source_density < ps_param.number_density && load_source_density < load_source_density_pre )
                {
                    zoom_factor = (1.0 - 0.5*itime/params->timestep_simulation) * (load_source_density_pre - load_source_density) / ps_param.number_density;
                    load_dn *= (1.0 + zoom_factor);
                    //cout<<"load_dn111  "<<load_dn<<"  "<<zoom_factor<<endl;
                }
                else if(load_source_density > ps_param.number_density && load_source_density > load_source_density_pre)
                {
                    zoom_factor = (1.0 - 0.5*itime/params->timestep_simulation) * (load_source_density - load_source_density_pre) / ps_param.number_density;
                    load_dn_temp = load_dn;
                    load_dn *= (1.0 - zoom_factor);
                    //cout<<"load_dn2222  "<<load_dn<<"  "<<zoom_factor<<endl;
                }

                if(load_dn <= 0.0)
                {
                    load_dn = load_dn_temp;
			    }
                */

                p_value = -(load_source_density - ps_param.number_density) / ps_param.number_density;
                i_value += p_value;

                if(load_source_density >= ps_param.number_density)
                {
                    d_value = (load_source_density - load_source_density_pre) / ps_param.number_density;
                }
                else
                {
                    d_value = -(load_source_density - load_source_density_pre) / ps_param.number_density;
                }

                load_dn_temp = load_dn;

                load_dn = load_dn_pid * (1.0 + ps_param.kp * p_value + ps_param.kd * d_value + ps_param.ki * i_value);

                if(load_dn <= 0.0)
                {
                    load_dn = 0.0;
                }
            }


            load_source_density_pre = load_source_density;


            load_number_double = load_dn * load_step * params->dt / params->species_param[ ps_param.species_number ].weight;
            load_number_rem_tot += load_number_double;
            load_number = load_number_rem_tot;
            load_number_rem_tot -= load_number;
            if(load_number < load_number_init_double)
			{
                load_dq_cell = 1.5 * load_temperature_init * (load_number_init_double - load_number);
                load_temperature = load_temperature_init;
			}
			else
			{
                load_dq_cell = 0.0;
                load_temperature = load_temperature_init * load_number_init_double / load_number;
			}
            //cout<<"nq end"<<endl;

        }
    }


    //for the case: input dn is vector in time, like Elms source
    if(loadTimeStepVector.size() > 1)
    {
        if(nextTimeStep == 0)
        {
            nextTimeStep = loadTimeStepVector[1];
        }
        if(itime > nextTimeStep)
        {
            load_temperature = load_temperatureVector[nextTimeStep_index];
            load_dn = load_dnVector[nextTimeStep_index];
            if(nextTimeStep_index < loadTimeStepVector.size() -1 )
            {
                //cout<<"ParticleSource1D_Load: nextTimeStep_index "<<nextTimeStep_index<<endl;
                nextTimeStep_index++;
                nextTimeStep = loadTimeStepVector[nextTimeStep_index];
            }
            else
            {
                nextTimeStep *= 1000;
            }

            load_step = 1.0 + load_number * params->species_param[ps_param.species_number].weight / (load_dn * params->dt);
            load_number_rem = load_dn * load_step * params->dt / params->species_param[ps_param.species_number].weight - load_number;
        }

    }


    if(ps_param.load_type == "nT" && load_cell_end[0] != load_cell_start[0])
    {
        for(int i = load_cell_start[0]; i <= load_cell_end[0]; i++)
        {
            cell_length = params->cell_length;
            cell_min[0] = params->local_min[0] + i * params->cell_length[0];

            auto cell_particles = species1D->particles_[i];
            cell_particles->reinitialize_particles(params, species1D->species_param, cell_min, cell_length, ps_param.mean_velocity, load_temperature, number_particle_in_each_cell[i]);
        }

    }
    else if(ps_param.load_type == "dn" && itime%load_step == 0 && load_cell_end[0] != load_cell_start[0])
    {
        load_number_rem_tot += load_number_double;
        load_number = load_number_rem_tot;
        load_number_rem_tot -= load_number;

        for(int i = load_cell_start[0]; i <= load_cell_end[0]; i++)
        {
            cell_length = params->cell_length;
            cell_min[0] = params->local_min[0] + i * params->cell_length[0];

            auto cell_particles = species1D->particles_[i];
            cell_particles->add_particles(params, species1D->species_param, cell_min, cell_length, ps_param.mean_velocity, load_temperature, load_number);
        }
    }
    else if(ps_param.load_type == "nq" && itime%load_step == 0 && load_cell_end[0] != load_cell_start[0])
    {
        for(int i = load_cell_start[0]; i <= load_cell_end[0]; i++ )
        {
            cell_length = params->cell_length;
            cell_min[0] = params->local_min[0] + i * params->cell_length[0];

            auto cell_particles = species1D->particles_[i];
            cell_particles->add_particles(params, species1D->species_param, cell_min, cell_length, ps_param.mean_velocity, load_temperature, load_number);
            if(load_dq_cell > 0.0)
            {
                //cell_particles->heat(params, species1D->species_param, load_dq_cell / (1.5 * cell_particles->size()) );
                //cell_particles->heat_new(params, species1D->species_param, ps_param, load_dq_cell);
                cell_particles->heat_new2(params, species1D->species_param, ps_param, load_dq_cell, ps_param.heat_number_ratio);
            }

            log_info<<"[ParitcleSource1D_Load.cpp] laod number "<<load_number;
            log_info<<"[ParitcleSource1D_Load.cpp] heat number "<<cell_particles->size() * ps_param.heat_number_ratio;
            log_info<<"[ParitcleSource1D_Load.cpp] load energy "<<1.5 * load_number * load_temperature;
            log_info<<"[ParitcleSource1D_Load.cpp] heat energy "<<load_dq_cell;
        }
    }
}

